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Tropical rainforests are significant contributors to the global cycles of energy, water and carbon. As a result, 
monitoring of the vegetation status over regions such as Amazonia has been a long standing interest of Earth 
scientists trying to determine the effect of climate change and anthropogenic disturbance on the tropical eco- 
systems and its feedback on the Earth's climate. Satellite-based remote sensing is the only practical approach 
for observing the vegetation dynamics of regions like the Amazon over useful spatial and temporal scales, but 
recent years have seen much controversy over satellite-derived vegetation states in Amazonia, with studies 
predicting opposite feedbacks depending on data processing technique and interpretation. Recent results 
suggest that some of this uncertainty could stem from a lack of quality in atmospheric correction and cloud 
screening. In this paper, we assess these uncertainties by comparing the current standard surface reflectance 
products (MYD09, MYD09GA) and derived composites (MYD09A1, MCD43A4 and MYD13A2 — Vegetation 
Index) from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Aqua satellite to results 
obtained from the Multi-Angle Implementation of Atmospheric Correction (MAI AC) algorithm. MAIAC uses a 
new cloud screening technique, and novel aerosol retrieval and atmospheric correction procedures which are 
based on time-series and spatial analyses. Our results show considerable improvements of MAIAC processed 
surface reflectance compared to MYD09/MYD13 with noise levels reduced by a factor of up to 10. Uncertainties 
in the current MODIS surface reflectance product were mainly due to residual cloud and aerosol contamination 
which affected the Normalized Difference Vegetation Index (NDVI): During the wet season, with cloud cover 
ranging between 90% and 99%, conventionally processed NDVI was significantly depressed due to undetected 
clouds. A smaller reduction in NDVI due to increased aerosol levels was observed during the dry season, with 
an inverse dependence of NDVI on aerosol optical thickness (AOT). NDVI observations processed with MAIAC 
showed highly reproducible and stable inter-annual patterns with little or no dependence on cloud cover, and 
no significant dependence on AOT (p<0.05). In addition to a better detection of cloudy pixels, MAIAC obtained 
about 20-80% more cloud free pixels, depending on season, a considerable amount for land analysis given the 
very high cloud cover (75-99%) observed at any given time in the area. We conclude that a new generation of at- 
mospheric correction algorithms, such as MAIAC, can help to dramatically improve vegetation estimates over trop- 
ical rain forest, ultimately leading to reduced uncertainties in satellite-derived vegetation products globally. 

© 2012 Elsevier Inc. All rights reserved. 


1. Introduction 

The Amazon basin encompasses almost half the tropical rainforest 
of the planet, storing an estimated 100 billion tonnes of carbon 
(Atkinson et al., 2011) and accounting for about 15% of global photo- 
synthesis (Malhi et al., 2008), while hosting roughly a quarter of the 
world's terrestrial species (Dirzo & Raven, 2003). Evapotranspiration 
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over Amazonia is a major factor of global energy and water circulation 
and affects precipitation rates across the Americas and large parts of 
the northern hemisphere (Werth & Avissar, 2002). Alterations of 
this ecosystem due to anthropogenic disturbance and climate change 
could have dramatic effects on the global carbon and energy balance 
(Gedney & Valdes, 2000). For instance, increased sea surface temper- 
atures in the Pacific ocean could intensify El Nino southern oscillation 
events and associated periodic Amazon droughts, while higher Atlantic 
sea surface temperatures and the northwest displacement of the inter- 
tropical convergence zone (Li et al., 2006) could cause these droughts to 
occur more frequently (Malhi et al., 2008). Some global circulation 
models (GCMs) predict that these drought effects alone may turn the 
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currently estimated carbon sink of Amazonia into a source, thereby sub- 
stantially affecting rates of global climate change (Lewis et al., 2011 ). 

As a result, accurate monitoring of the Amazon is essential for es- 
timating future climate scenarios, predicting global energy fluxes and 
precipitation rates, and assessing the success of conservation efforts, 
such as the United Nation's REDD (reducing emissions from defores- 
tation and degradation) initiative (Mollicone et al., 2007). While sat- 
ellite remote sensing provides the only viable means to observe the 
Amazon in a spatially comprehensive and temporally frequent fash- 
ion, spaceborne optical assessment of tropical vegetation is inherent- 
ly difficult due to high cloud cover, high aerosol concentration from 
biomass burning and limited physical access to ground validation 
sites. Consequently, recent years have seen much controversy over 
satellite derived estimates of the vegetation status in Amazonia. For 
instance, Saleska et al. (2007) reported an increase in vegetation 
greenness over the Amazon during the 2005 drought based on 
MYD13A2 observations from MODIS. This study was subsequently 
challenged by others, most notably Samanta et al. (2010) and 
Atkinson et al. (201 1 ) on the basis of poor data quality and processing 
methodology. In response to a second drought event in 2011, Xu et al. 
(2011) reported a widespread decline in Amazon greenness, thereby 
directly contradicting the results reported by Saleska et al. (2007) for 
the 2005 drought. Similarly, Zhao and Running (2010) reported a 1% 
decrease in global primary productivity between 2001 and 2009, 
which was based largely on a decrease in primary productivity (and 
fpAR ) in the Amazon, a result subsequently challenged by Samanta et 
al. (2011) and Medlyn (2011). 

Recent work (Forster & Ramaswamy, 2007; Levy et al., 2010) sug- 
gests a substantial uncertainty of atmospheric aerosol properties from 
MODIS data over the Amazon region as potential cause of these discrep- 
ancies (Zelazowski et al., 2011). These uncertainties could introduce 
errors in discrete and continuous land surface parameters, such as 
land cover classes (Lu & Weng, 2007) and vegetation indices (Myneni 
& Asrar, 1994) throughout the Amazon basin (Zelazowski et al., 2011). 
The current MODIS land surface reflectance product (MOD09/MYD09) 
routinely corrects for atmospheric scattering using models of the radia- 
tive transfer of light through the atmosphere (Vermote & Kotchenova, 
2008). While this technique is widely applied, it uses a single observa- 
tion to characterize a pixel value driven by two main unknowns, aerosol 
optical thickness (AOT) and surface reflectance (SR) and, as a result, a 
priori assumptions are required to describe their relationship. For in- 
stance, the current aerosol retrieval algorithm of MOD04/MYD04 
(Kaufman et al., 1997; Levy et al., 2007) relates surface reflectance in 
the visible (blue and red) spectral bands with MODIS band 7 (2.1 pm) 
reflectance by using a prescribed spectral regression coefficient (SRC). 
A Lambertian surface model is then applied for aerosol retrievals and 
atmospheric correction. While this technique greatly simplifies data 
processing, the Lambertian assumption reduces the anisotropy of the 
derived reflectance and introduces an error that depends on the aerosol 
amount and the view-observer geometry (Lyapustin, 1999; Wang et al., 
2010 ). 

In addition to the uncertainties from the Lambertian assumption, 
pixel-based algorithms have to rely on spectral reflectance and thermal 
thresholds for cloud masking. However, lack of a priori knowledge 
about the specific cloud free reflectance makes the distinction between 
cloudy and clear observations difficult (Yang & Di Girolamo, 2008). Tra- 
ditionally, a generic land type classification has been used to address 
this limitation by substituting a standard reflectance value for clear 
pixels; however, considerable uncertainties remain due to the wide 
natural variability of both land surface and cloud reflectance (Lyapustin 
et al., 2008; Rossow & Garder, 1993). 

MAIAC grids MODIS LIB data to a 1 km resolution, and accumu- 
lates measurements of the same surface area from different orbits 
(view geometries) for up to 16 days of observations for equatorial 
and up to 4 days for polar regions using a moving window approach. 
The MAIAC cloud mask (CM) algorithm composes a dynamically 


updated reference clear-sky image of the surface from spatial and 
time series analyses. The knowledge of reference clear-sky reflectance 
in addition to spectral and thermal reflectance tests (Ackerman et al., 
1998) has been shown to improve cloud detection (for detail, see 
Lyapustin et al., 2008). MAIAC aerosol retrieval (Lyapustin et al., 
2011b) and atmospheric correction (Lyapustin et al., this issue) algo- 
rithms use an advanced radiative transfer theory with anisotropic 
land surface reflectance (Lyapustin & Knyazikhin, 2001 ) parameterized 
by the Ross-Thick Li-Sparse (RTLS, Roujean et al., 1992) bidirectional 
reflectance model (Lyapustin et al., 2011a). Once RTLS in MODIS band 
7 (2.13 pm) is initialized, the aerosol algorithm derives a spectral re- 
gression coefficient (SRC) for each 1 km grid cell by simultaneously 
processing 25x25 km 2 blocks of pixels from four or more cloud-free 
observations obtained from different view angles. Knowledge of SRC 
allows aerosol retrievals at 1 km resolution over spatially heteroge- 
neous surfaces, and helps to avoid dependence of derived AOT on sur- 
face brightness (Lyapustin et al., 2011b). Finally, the atmospheric 
correction algorithm (see Lyapustin et al., this issue) derives parameters 
of the RTLS model for each 1 km grid cell along with spectral surface al- 
bedo and a bidirectional reflectance factor (BRF). An overview of all 
components of MAIAC is provided in Lyapustin and Wang (2009). 

Algorithms like MAIAC, with more complete descriptions of the 
physical system and limited role of empirical assumptions, hold 
promise to overcome many of the restrictions currently faced by con- 
ventional satellite retrievals over tropical regions. In this paper, we 
compare MODIS standard atmospheric correction from Aqua (using 
the MYD09 and MYD13 reflectance products) to MAIAC results over 
the Amazon basin and assess the potentials and limitations of satellite 
retrievals of vegetation greenness over Amazonia. The implications 
for determining the biophysical state of the rainforest and detecting 
changes over time are being discussed for the different products. 

2. Methods 

2.1. Data 

The study encompasses two MODIS tiles (hi 1 v09 and hi 1 v08), an 
area of 2.88 million km 2 , spanning 10°N to 10°S in latitude and 70°30' 
W to 62°W in longitude (Fig. 1 ). The area is characterized by seasonal 
tropical savannah and seasonal rain forest in the north, while the 
southern part consists of evergreen tropical forest. We used MODIS 
data from the Aqua satellite platform collected between 2002/07/04 
and 2010/12/31. Collection 5 data of the MYD09GA 1 km daily surface 
reflectance product and the MYD09A1 eight-day composite product 
were obtained from EOS data gateway of NASA's Goddard Space 
Flight Center (https://wist-ops.echo.nasa.gov/api/) and data mosaick- 
ing was performed using the MODIS reprojection tool (MRT). Clouds 
were masked by means of the 'state_lkm' scientific dataset (SDS) in- 
cluded in the MYD09GA product ('state_500m' in case of MYD09A1), 
which is based on two cloud detection algorithms, the MOD/MYD35 
cloud mask (Frey et al., 2008) and an additional, internal cloud screen- 
ing (Vermote et al., 2008). In addition to cloud masking, all MYD09 data 
were quality filtered using the MYD09 quality (Qa) flags and only cloud 
free pixels with high data quality were passed and used for all subse- 
quent analyses (Vermote et al., 2008). An overview of the quality and 
cloud flags set for quality assurance is given in Table 1. 

In addition to MYD09, we also acquired the MODIS 16-day 500 m 
vegetation index (VI) product, MYD13A2, for the same observation pe- 
riod from the same data gateway. Data quality was assured using the 
"500 m 16 days VI Quality’ flags as well as the ‘500 m 16 days pixel re- 
liability 4 layer to allow only pixels with good quality and high reliability 
rating to pass (see Table 1 ). The MYD13A2 product supplies both NDVI 
and the Enhanced Vegetation Index (EVI), however in this study we 
focussed on the effect of atmospheric correction and cloud screening 
on NDVI. The compositing algorithms, MYD09A1 and MYD13A2, accu- 
mulate observations over 8 and 16 days, respectively, and use pixel 


372 


T. Hilker et al. / Remote Sensing of Environment 127 (2012) 370-384 


V 

Y 


70°0 , 0"W 67°30 , 0"W 65°0 , 0 ,, W 62°30 , 0"W 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 
17 



NDVI 



70 o 0'0"W 67°30'0 ,, W 


65°0'0"W 62°30'0"W 

0 75 150 300 450 600 

■ ■ ■ Kilometers 


Fig. 1. Study site. The study area encompasses two MODIS tiles, spanning an area of 2.88 million km 2 . 


values with darkest surface reflectance (highest NDVI in case of 
MYD13A2) and closest to nadir view angle to represent a given lo- 
cation and compositing cycle. 


MODIS Aqua calibrated and geometrically corrected (Level IB) 
data acquired between 2002/07/04 and 2010/12/31 for the same 
area were also processed by MAIAC. For this study, we used the 


Table 1 

Quality /cloud flags of the MYD09GA, MYD09A1, MYD13A2 and MCD43A2 (for MCD43A4 quality) products used in this study. Flags that are not mentioned were not used. 


Product 

SDS name 

Flag 

Accepted values 


MYD09GA 

QC 500 m 

MODLAND QA bits 

00 

(ideal quality — all bands) 



band 1 data quality 

0000 

(highest quality) 



band 2 data quality 

0000 

(highest quality) 



atm. corr. performed 

1 

(yes) 


State 1 km 

cloud state 

00 

(clear) 



cloud shadow 

0 

(no) 



aerosol quantity 

00/01 

(climatology/low) 



cirrus detected 

00 

(none) 



internal cloud flag 

0 

(no cloud) 



fire flag 

0 

(no fire) 



pixel adjacent to cloud 

0 

(no) 

MYD09A1 

QC 500 m 

see ‘QC 500 m’ for MYD09GA 




State 500 m 

see ‘State 1 ’ km for MYD09GA 



MYD13A2 

VI Quality 

MODLAND QA Bits 

00 

(VI produced with good quality) 



VI Usefulness 

0000 

(Highest quality) 




0001/0010/0100/1000 

(Lower quality) 



aerosol quantity 

00/01 

(climatology/low) 



adjacent cloud 

0 

(no) 



mixed cloud 

0 

(no) 



possible shadow 

0 

(no) 


Pixel reliability 

Rank Key 

0 

(Good data — use with confidence) 

MCD43A2 ( appendix ) 

Albedo Quality 


0 

(Processed, good quality) 




1 

(Processed, see other QA) 


Albedo Band Quality 

Band 1 

0000 

(best quality) 




0001 

(good quality) 



Band 2 

0000 

(best quality) 




0001 

(good quality) 



T. Hilker et al. / Remote Sensing of Environment 127 (2012) 370-384 


373 


following MAIAC outputs: cloud mask, aerosol optical thickness 
(AOT), bidirectional reflectance factors (BRF, also known as surface 
reflectance), and RTLS model parameters for the red and near-IR 
MODIS bands. MAIAC data does not include quality flags, but rather 
provides observations only for those pixels that were found not to 
be cloud contaminated and for which a stable aerosol correction 
could be made. 

We used the BRF (or surface reflectance) values from both the 
MYD standard products and MAIAC to compare the effects of cloud 
masking and aerosol correction on the assessment of vegetation sta- 
tus from MODIS. A separate assessment of the BRDF-normalization 
(to remove view geometry variations) on MODIS NDVI can be found 
in the appendix. 

2.2. Approach 

We conducted a comparative analysis of cloud mask, aerosol, 
and NDVI between MAIAC and the respective operational MODIS 
algorithms. Performance of the two cloud screening algorithms 
(MYD35/MYD09 and MAIAC) was assessed relative to each other, by 
computing errors of omission/commission of MYD35/MYD09 using 
MAIAC as a reference. Further, MAIAC aerosol retrievals were qualita- 
tively compared to the MODIS operational aerosol product MYD04 as 
the internal aerosol product of MOD09, which is based on similar prin- 
ciples as the Dark Target method of MYD04 (Vermote & Kotchenova, 
2008) but at a 1 km resolution, is not reported. The MYD04 Level 2 
data is available at a 10 km spatial resolution and was gridded to the 
same sinusoidal projection as MAIAC, MYD09 and MYD13 to allow a 
direct comparison between the data products. 

The normalized difference vegetation index was computed using 
MODIS bands 1 and 2 from the two MYD09 products (MYD09GA 
and MYD09A1 ) and from MAIAC, and was also acquired directly 
from MYD13A2. The MAIAC NDVI was computed from MAIAC BRF 
data in order to be directly comparable to MYD09GA, MYD09A1 and 
MYD13A2 NDVI. The area studied in this work includes seasonal trop- 
ical forests and savannahs in the north (Fig. 1 ), which are expected to 
show seasonal variability in vegetation greenness driven by the regu- 
lar changes between dry and wet season. Relatively little variability is, 
however, expected in the southern part of the study area, which is 
covered by dense evergreen forest. As a result of the climatological 
conditions found in the study area, it can be assumed that vegetation 
indices related to greenness should be relatively stable over short pe- 
riods of time (a few days) as there is no physiological reason for rapid 
and notable changes vegetation greenness (>5-10%). The short-term 
stability of NDVI can therefore be seen as an indicator of data quality, 
since high frequency changes in NDVI will most likely be driven by 
cloud and aerosol artifacts rather than changes in vegetation status. 

3. Results 

Fig. 2 shows a comparison of NDVI time series obtained from 
MYD09GA, MYD09A1, MYD13A2 and MAIAC across six 50x50 km 2 
spatial subsets of the study area. The subsets were selected using reg- 
ularly spaced raster sampling to obtain representative regions for 
both seasonal and non-seasonal vegetation (as marked in the map). 
To allow a comparison between the daily, biweekly and monthly 
products, we first computed monthly mean values for each pixel by 
averaging all valid (that is cloud screened, and quality filtered) pixels 
acquired during a month. The arithmetic means and standard devia- 
tions shown in Fig. 2 were then obtained by averaging the monthly 
values across each 50x50 km 2 block. It is important to note that our 
processing approach excludes temporal variability in the standard devi- 
ation shown since each pixel is represented by only one monthly value, 
regardless of whether MYD09GA, MYD09A1, MYD13A2 or MAIAC was 
used. As a result, the standard deviations describe only variability due 
to heterogeneity of the landscape. Differences between these standard 


deviations can be attributed to differences in cloud screening and atmo- 
spheric correction since all datasets describe the identical area over the 
same time period. The largest variability in NDVI was found for the 
MYD09GA 1 km daily surface reflectance product across all samples, 
while the errorbars were smaller for the MYD09A1 composite (Fig. 2). 
The 16-day composite showed the least variability of all conventionally 
processed products, but scattering was still about twice as large as that 
from MAIAC dataset, which showed an up to 1 0-fold reduction in spatial 
variability compared to the daily MYD product. The northern sites 
exposed a clear seasonal signal of NDVI across all examined MODIS 
datasets, while little seasonal variability was found for the two equato- 
rial and the southern regions. 

3.1. Cloud masking 

Overall, a good correlation was found between MYD35 and MAIAC 
in cloud detection (r 2 = 0.75), however, MAIAC cloud screening 
yielded significantly more cloud free observations than MYD35/ 
MYD09 (Figs. 3, 4). Fig. 3A presents a false color NIR image of the en- 
tire study area observed on July 16, 2002 (day of year, DOY=197); 
the high cloud fraction typical for the Amazon basin is visible in 
both cloud screening products (Fig. 3B and C). For this overpass, the 
MYD09/MYD35 cloud mask (Fig. 3B) identified 6.84% of the image 
as cloud free, significantly less than the MAIAC algorithm (Fig. 3C), 
which found cloud free pixels in 14.34% of the entire image area. 
Depending on algorithm and season, the detected cloud cover was 
between 75% and 99% of the total area across all examined MODIS 
scenes (Fig. 4). Seasonal variation in cloud cover was in the order of 
10-15% between the wet and dry seasons. On average, the area 
marked as cloudy by MAIAC was about 4% less than that masked by 
MYD09/MYD35 which, given the cloud cover, provides on average 
between 20-80% more cloud-free data for the land analysis in Amazon 
basin, depending on season. 

Fig. 5 shows a comparative metric of cloudy /clear pixel classifica- 
tion of MYD09/MYD35 relative to that of MAIAC. The red line shows 
the relative accuracy of the MYD09/MYD35 cloud detection (the 
ratio of the number of detected cloud pixels to that of MAIAC is called 
user accuracy). The difference between 100% and this line represents 
the error of commission (relative to MAIAC), when MYD09/MYD35 
classified a pixel as being cloudy, whereas the reference algorithm 
(MAIAC) did not. As the majority of pixels in each scene is classified 
as cloudy in both algorithms, the relative error of commission is low 
(~6%). The difference between 100% and the black line represents 
the error of omission when MYD09/MYD35 failed to detect a cloud, 
relative to MAIAC cloud mask. 

The relative accuracies for cloud-free pixels derived from MYD09/ 
MYD35 and MAIAC were significantly lower, as only 1% to 25% of each 
image area was cloud-free, and as a result, a misclassification error of 
the same number of pixels is relatively larger. The green and blue lines 
show the omission/commission error for the clear pixels. The green line 
gives the percentage of pixels classified as clear by MYD09/MYD35 
while MAIAC classified them as cloudy, whereas the blue line shows 
cases where MYD09/MYD35 classified a pixel as cloudy while MAIAC 
did not. 

3.2. Aerosol optical depth 

To illustrate differences in aerosol optical thickness between the 
MAIAC and operational MODIS aerosol retrievals (MYD04), we 
contrasted the distribution of AOT for a seasonal (northern) and a 
non-seasonal (southern) region of the study area using the year 
2005 as an example; a year that saw a significant drought with asso- 
ciated intense biomass burning (Fig. 6). The selected areas corre- 
spond to the top left (Fig. 6A and C) and bottom left regions (Fig. 6B 
and D) shown in Fig. 2. To account for difference in resolution 
(1 km in MAIAC vs 10 km in MYD04), data were averaged over a 
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Fig. 2. Monthly observations of NDVI from MYD09GA (1 km daily surface reflectance), MYD09A1 (8-day composites), MYD13A2 (16-day vegetation index product) and MAIAC 
obtained over a 50 x 50 km area as indicated in the map. The errorbars represent the mean spatial variability over the 50 x 50fkm subset (Please note that the points and errorbars 
of the products are slightly offset from each other for purposes of readability). 


larger area of 500 x 500 km 2 . The top row in Fig. 6 illustrates MAIAC 
aerosol optical thickness at 0.47 pm, the bottom row shows the 
corresponding MYD04 AOT at 0.55 pm. In the northern region (left 
column), AOT was relatively stable throughout the year with most 
observations around 0.2-0.3 and much less frequent higher AOT 
values. In the southern region (right column), AOT showed a clear 
seasonality with low values during the wet season and high average 
AOT related to biomass burning during the dry season. A somewhat 
lower MYD04 AOT is mostly explained by the wavelength difference 
of the AOT product (Levy et al., 2007; Lyapustin & Wang, 2009). 
While MAIAC and MYD04 exposed similar patterns for both regions, 
the MAIAC AOT product shows a smoother frequency distribution 
over time (y-axis), in part due to higher spatial resolution. Fig. 7 
shows a direct comparison between AOT derived from MAIAC and 
MYD04 across the entire study area (the two tiles), resampled to 
100 km 2 to account for the difference in spatial resolution existing 
between the two products. AOT derived from MAIAC shows a good 
correspondence to that derived from MYD04 throughout the study 
period (data shown between 2002 and 2010). The relationship was 
not significantly biased with an r 2 of 0.84 at p<0.01. 

3.3. NDVI 

To analyze the implications of the observed differences in cloud 
masking and AOT on vegetation studies, we assessed the spatial and 
temporal variability in NDVI using the same northern and southern 


sample areas shown in Fig. 2. Figs. 8 and 9 present NDVI estimates 
(along z-axis) as a function of time during 2002-2010 (x-axis). 
The y-axis shows the spatial variability of NDVI when averaged over 
different areas, from 2x2 to 50x50 km 2 . The color represents the 
standard deviation (ct) of the spatially averaged NDVI estimates. 
Since the tropical landscape is fairly homogeneous (particularly in 
the south), a mostly characterizes processing errors from clouds 
and aerosols, at least over smaller areas (up to -100 km 2 ). The color 
bar shows that at these scales the MAIAC NDVI uncertainty is small, 
about 0.01-0.02, which represents the total retrieval error (Figs. 8A 
and 9A). The uncertainty increases to about 0.04-0.05 at larger scales 
(^30-50 km 2 ) reflecting the landscape spatial variability. The vege- 
tation signal in the north was strongly driven by seasonal effects 
(Fig. 8A), while in the south, almost no temporal or spatial changes 
were observed (Fig. 9A). Despite the differences in AOT, cloud cover 
and vegetation dynamics found between the northern and southern 
study area, MAIAC produced a highly consistent and reproducible 
NDVI time record. 

Compared to MAIAC, a of MYD09GA-based NDVI estimates 
was about 10-fold higher and ranged between 0.15 and 0.20 (>0.25 
in several areas) even when aggregating relatively few pixels 
(^100 km 2 ) (Figs. 8B and 9B). Results for the southern site show 
that the uncertainty in NDVI remained high at all spatial scales, 
indicating that the algorithm errors are significantly higher than 
spatial variability of the vegetation. This increased uncertainty in 
the MYD09 estimates compared to MAIAC can be attributed to both 
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Fig. 3. Comparison of MYD09/MYD35 and MAIAC cloud mask. (A) False color infrared image acquired over the study area on July 16, 2002. (B) MYD09/MYD35 derived cloud mask 
overlaid on top of the same image (C) MAIAC derived cloud mask. 


cloud masking and aerosol correction errors. The low AOT found at 
the northern site and also at the southern site during the wet season, 
at least in 2005, suggests that this variability is mainly due to the 
cloud masking errors. To assess the contributions of AOT and cloud 
masking in more detail, we repeated the same analysis as shown in 
Figs. 8B and 9B but using MYD09GA product filtered by the MAIAC 
cloud mask. The results of these analyses are shown in Figs. 8C 
and 9C, for the northern and southern study area, respectively. The 
resulting NDVI values were notably smoother and had a much more 
consistent time record approaching the MAIAC results (shown in 


Figs. 8A and 9A). This confirms that most of the MYD09GA errors 
were due to cloud masking. Compared to Figs. 8B and 9B, a was re- 
duced to 0.02 and 0.1, respectively, suggesting that cloud contamina- 
tion was responsible for about 80% of uncertainties seen in these 
Figures. Nonetheless, even with the MAIAC cloud filter applied, the 
MYD09GA NDVI uncertainty was about twice that of MAIAC. 

The MYD09A1 and MYD13A2 compositing (Figs. 8D and 9D, and 
Figs. 8E and 9E, respectively) had only limited effects on improving 
the uncertainties of the conventional NDVI product. For the northern 
study area, a ranged between 0.03 to 0.1 for the MYD09A1 product, 
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Fig. 4. Comparison between the fraction of cloud free pixels acquired throughout the study period from MAIAC (black) and MYD09/MYD35 (red). MYD09/MYD35 was consistently 
more conservative than MAIAC, cloud cover was heavily dependent on season. 
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Fig. 5. Errors of omission and commission for cloud and non cloud pixels acquired from confusion matrices of each acquisition date comparing MYD09/MYD35 cloud mask to 
MAIAC. One minus the red line represents the fraction of cases where MYD09/MYD35 classified pixel as cloudy, whereas MAIAC found this pixel to be cloud free. One minus the 
black line represents the fraction of cases, in which MAIAC identified a pixel as being cloudy, while MYD09/MYD35 marked it as cloud free. The green line shows the percentage 
of cloud free pixels found by MYD09/MYD35 but marked as cloudy in MAIAC, whereas the blue line shows cases in which MAIAC had identified a pixel as cloud free but not 
MYD09/MYD35. 


while for the 16-day MYD13A2 product, a was reduced to 0.025 and 
0.09 compared to MYD09GA. Similarly, for the southern area, a 
ranged between 0.02 and 0.1 for MYD09A1, and between 0.01 and 
0.05 for MYD13A2 products. Residual cloud contamination was ap- 
parent in the composite products, especially in the south, showing 


inconsistent annual patterns and frequent spatial outliers (compare 
also Fig. 2). 

The results presented in Figs. 8 and 9 are consistent with the find- 
ings presented in Fig. 10, which shows a direct comparison of NDVI 
obtained from MAIAC and MYD09GA across the entire study area. 



AOT AOT 


Fig. 6. Variation in aerosol optical thickness (AOT) over a northern (7°30’ N, 70° W , Figure A) and a southern (7°30' S, 70° W , Figure B) subset of the study area (size 500 x 500 km 2 
each). The x-axis represents the time (year 2005), they-axis shows AOT bins (width 0.01 ) and the z-axis represents counts per AOT bin (as fraction of total number of observation at 
a particular time). The color corresponds to the total number of observation at a particular time. Figures C and D show the same AOT estimates, but derived from MYD04 at 10 km 
resolution. 
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Fig. 7. Comparison between MAIAC AOT and MYD04 for the entire study area (MAIAC 
data were resampled to 10 km 2 to match the MYD04 resolution) between 2002 and 
2010, filtered by both, the MYD35 and MAIAC cloud masks. 


Shown NDVI values were accumulated over an eight day period of 
DOY 356-364, 2010 to collect more observations, and were placed 
on a regular grid using 3D spline interpolation. While the MAIAC 
(Fig. 9A) and MYD09GA (Fig. 9B) processed NDVI data both captured 
the phenological differences between the seasonal (north) and the 
non-seasonal (south) forest, a much higher spatial variability of 
NDVI was observed from MYD09GA (Fig. 10B) compared to MAIAC 
(Fig. 10A). This spatial variability is not supported by changes in the 
land cover type (Fig. 2) or heterogeneity in tropical vegetation and 
is therefore indicative of the noise level discussed above. While this 
Figure shows only one 8-day period, we created a movie file covering 
the entire study area for a two-year period. This movie, which is 
available in the supportive online material, provides a visual demon- 
stration of algorithm differences in assessment of the temporal vari- 
ability of NDVI. 

The increased noise levels limited the ability of MYD09 observa- 
tions to describe shifts in vegetation greenness. This can be demon- 
strated when comparing the distributions of NDVI observations 
across a latitudinal gradient from 10°N to 10°S by binning NDVI 
values into strata (based on a bin width of 0.01 ) and counting the per- 
centage of observations per bin and latitude (i.e. aggregating across 
all longitudes). For the non-seasonal tropical forest, the vast majority 
of NDVI observations are expected to be at the higher end of the scale 
(0.8-0.9) with little variability likely across southern latitudes, 
whereas seasonal vegetation in the north should show clear shifts in 
NDVI distributions when comparing dry and wet seasons. Fig. 1 1 pre- 
sents NDVI distributions for one day in the northern tropics wet 
season (Figure A and B) and a northern tropics dry season (Figure C 
and D). Rather than absolute counts of observations, percentages of 
total cloud free observations are shown along the z-axis in order to 
normalize for latitudinal variations in cloud cover (as cloud cover 
drives the number of clear pixels that can be aggregated across each 
longitude). MAIAC observed NDVI distributions where highly consis- 
tent across all southern latitudes with values ranging between 0.7 
and 0.9 (left column). A slightly lower range of NDVI values was 
found in the northern latitudes, due to lower vegetation cover in 
the tropical savannah regions (Fig. 11 A, compare Fig. 1). In the 
north, a clear reduction in NDVI values was observed during the dry 
season (Fig. 1 1C) while distributions in the south showed little or 
no change. The distribution of MYD09GA data (Fig. 11B and D) was 
in both cases more skewed towards lower NDVI values and the 


relative number of high NDVI values was lower and less equally dis- 
tributed across latitudes. 

MYD09 processed NDVI observations were also mostly distributed 
around 0.7 and 0.9 across southern latitudes, however, the distribu- 
tion was notably wider and skewed towards the lower range (<0.7), 
likely because NDVI observations were depressed by the cloud and 
aerosol contamination. In addition, the NDVI distribution varied 
strongly even across southern latitudes, which is not supported by 
latitudinal changes in landcover and can therefore be attributed to 
variations in cloud and aerosol contamination. Similarly, the shifts 
between dry and wet season the northern latitudes are much less ap- 
parent in MYD09 compared to MAIAC, thus showing a more limited 
ability of MYD09 to detect seasonal changes in greenness in condi- 
tions of high cloudiness. 

Finally, we also analyzed the dependence of NDVI on AOT directly, 
by using the two 10x10 km 2 sample areas earlier presented in Figs. 7 
and 8. In order to eliminate potential differences due to cloud 
masking and quality screening and to single out the aerosol effect, 
data shown in Fig. 12 were filtered using both MAIAC and MYD09/ 
MYD35 cloud mask and quality flags. Data points acquired between 
2002 and 2010 are shown; the color of the data points corresponds 
to the respective month of acquisition. MYD09 processed data 
showed a statistically significant relationship with both MAIAC and 
MYD04 derived AOT (r 2 = 0.2 p<0.05, Fig. 12B and D and r 2 = 0.15 
p<0.05, Fig. 12C and E), as high AOT values depressed NDVI observa- 
tions in both northern and southern sites. The random distribution of 
the color codes shows that this dependence cannot be explained by 
seasonality of NDVI and AOT, suggesting a relatively high aerosol 
effect on MYD09 NDVI. MAIAC NDVI showed only a weak AOT depen- 
dence in both regions, which was not, however, statistically signifi- 
cant at p<0.05 (Figures A and D). 

4. Discussion 

This study presented a comparative analysis between the current- 
ly operational and a new atmospheric correction algorithm MAIAC to 
monitor the vegetation status in tropical regions from MODIS data. 
We have shown that cloud screening and correction for aerosols are 
critical for reducing uncertainties of remote sensing of vegetation 
over tropical regions which confirms results obtained by other groups 
(Di Rosa et al., 2009). For instance, Figs. 8 and 9 suggest that cloud 
leakage was the single most important factor driving high spatial 
and temporal variability in MYD09GA derived NDVI. The time series 
approach implemented in MAIAC helped to significantly reduce 
these uncertainties as it was able to detect clouds more reliably 
while providing a significantly larger amount of cloud free observa- 
tions. This resulted in MAIAC showing stable and well-reproducible 
patterns in vegetation greenness with errors and uncertainties re- 
duced by a factor of up to ten compared to current state-of-the-art 
retrievals. 

While we did not conduct a detailed analysis here, the main im- 
provement in MAIAC cloud screening is probably achieved for low 
and thin or small (sub-pixel) clouds, which are partly captured by 
the MAIAC aerosol retrievals (Lyapustin et al., 2011b, 2012) but oth- 
erwise are hard to detect using conventional spectral or brightness 
temperature tests. It has been a common land community practice 
to mitigate the effect of missed clouds by compositing data over 
8-day (MYD09A1) or 16-day periods (MYD13A2). This effectively im- 
plements an additional cloud filter, as values with lower reflectance 
are preferred over higher reflectance values, thereby eliminating 
data with a higher potential for cloud contamination. Our results 
presented in Figs. 2, 8 and 9, show that while this technique helped 
to reduce ct somewhat, the composited data still showed high levels 
of uncertainties over the study area, making analysis of vegetation 
greenness or trend predictions over Amazonia challenging. For in- 
stance, obtaining accurate results from conventional MODIS products 
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Fig. 8. Comparison between MAIAC (A) and MYD09 (B) derived NDVI values aggregated over a range of different area sizes from 4 to 2500 km 2 for the northern subset of the study 
area (7°30’ N, 70° W). The x-axis shows the size of the aggregated area (pixel), the y-axis represents the time (2002-2010) the z-axis shows the corresponding NDVI value. The 
color of the surface corresponds to the standard deviation obtained from averaging over different areas. Figure C shows MYD09 derived data, but using the MAIAC cloud mask. 
Figure D shows the same product derived from the 8 day MYD09A1 composite, Figure E represents the MYD13A2 16 day product. Please note that the number of observations for 
these composites is less compared to the daily products. 


(MYD09GA, MYD09A1 , MYD13) requires large statistics and very 
skillful and sophisticated data filtering approach (e.g., Samanta et 
al, 2012). 

While additional data quality analysis and filtering may improve 
the results, it also substantially reduces the amount of information. 
For instance, composited data yielded better results than the daily 
MODIS products (Figs. 2, 8, 9, Al, A2) but less than half as many pixels 


were observed over a 16 day period than observed by MAIAC (see for 
instance Fig. A3) over the same period of time, not including the 
higher number of MAIAC observations for the same pixels during 
the compositing period. Compared to the standard product outputs, 
the MAIAC derived cloud mask yielded on average about 25% more 
cloud free observations, which is a significant gain considering that 
between 75% and 99% of the images was covered by clouds at any 
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Fig. 9. Comparison between MAIAC (A) and MYD09 (B) derived NDVI values aggregated over a range of different area sizes from 4 to 2500 km 2 for the southern subset of the study 
area (7°30 / S, 70° W). The x-axis shows the size of the aggregated area (pixel), the y-axis represents the time (2002-2010), the z-axis shows the corresponding NDVI value. The 
color of the surface corresponds to the standard deviation obtained from averaging over different areas. Figure C shows MYD09 derived data, but using the MAIAC cloud mask. 
Figure D shows the same product derived from the 8 day MYD09A1 composite, Figure E represents the MYD13A2 16 day product. 


given time. The potential for increasing the number of cloud free ob- 
servations is an important finding, as it can help to enhance the fre- 
quency of satellite based assessment of forest loss and landscape 
disturbance over tropical regions (Hansen et al., 2008, 2009) and 
allow more timely and accurate predictions of the status of the 
remaining tropical forest. 


Consistent with previous work, our results have shown that atmo- 
spheric effects can depress NDVI observations (Myneni & Asrar, 1994; 
Xiao et al., 2003). The relationship between AOT and MYD09 NDVI 
presented in Fig. 12 suggest that MYD09GA collection 5 surface re- 
flectance contained residual atmospheric effects as both northern 
and southern study sites showed a noticeable AOT-dependence of 
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Fig. 10. Spatial distribution of NDVI over the entire study area using NDVI observations processed with MAIAC (A) and MYD09 (B) acquired over a 10 day period between DOY 355 
and DOY 364, 2010. The x- and y-axis represent geographical longitude and latitude, respectively, the z-axis shows NDVI. A three dimensional spline function has been fitted 
through all valid NDVI observation acquired by either algorithm over the 10 day period to obtain a continuous surface from the discrete observations. 




Fig. 11. Distribution of NDVI across a latitudinal gradient from 10°S to 10°S for a day in late northern latitude summer (top row) and late northern latitude winter (bottom row). The 
x-axis represents NDVI bins (bin size = 0.01), the y-axis shows the latitudinal range. The z-axis represents NDVI observations per bin (aggregated across all longitudes, data are 
shown as fraction of the total number of observations at a particular latitude). The left column (Figure A and C) shows MAIAC derived data, the right column (Figures B and D) 
shows MYD09 derived observations. The total number of observations per bin is indicated by the color. 
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Fig. 12. Relationship between NDVI and MAIAC derived aerosol optical thickness (AOT) at a northern (7°30' N, 70° W) and a southern (7°30' S, 70° W) subset of the study area (size 
100 km 2 each). The left column (Figure A and D) shows MAIAC processed observations, the central column (Figure B and E) shows NDVI obtained from MYD09 and AOT derived 
from MAIAC. The right column (Figure C and F) shows MYD09 derived NDVI compared to MYD04 derived AOT. Data acquired over the full 8 year period are shown (2002-2010). 
The marker color corresponds to the month of acquisition as indicated by the bars on the right. 
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Fig. Al. Monthly observations of NDVI from MCD43A4 (1 km 16-day NBAR view angle normalized), and MAIAC BRFn obtained over a 50x50 km area as indicated in the map. The 
errorbars represent the mean spatial variability over the 50 x 50 km subset (Please note that the points and errorbars of the products are slightly offset from each other for purposes 
of readability). 
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Fig. A2. Comparison between MAIAC BRFn (left column A, C) and MCD43A4 (right column B,D) derived NDVI values aggregated over a range of different area sizes from 4 to 
2500 km 2 for the northern (7°30 / N, 70° W, top row) and southern subset of the study area (7°30’ S, 70° W, bottom row). The x-axis shows the size of the aggregated area 
(pixel), the y-axis represents the time (2002-2010) the z-axis shows the corresponding NDVI value. The color of the surface corresponds to the standard deviation obtained 
from averaging over different areas. 


MYD09 NDVI. Partly, it can be explained by the Lambertian surface as- 
sumption used in MYD09, which does not discriminate between the 
direct and diffuse irradiance, the effective surface reflectance for 
which is different. The relative amount of direct vs. diffuse radiative 
flux at the surface depends on the aerosol amount and path through 
the atmosphere (solar and view zenith angle). This causes biases in 
the surface reflectance which are dependent on the viewing geometry 


and AOT even with perfect knowledge of AOT (Lyapustin, 1999; Wang 
et al., 2010). 

This study has shown that MAIAC algorithm can yield dramatically 
enhanced estimates of commonly derived vegetation indices such as 
NDVI over tropical regions. As a result, MAIAC may allow new oppor- 
tunities for re-assessing the state of terrestrial ecosystems even in 
very cloudy regions like the Amazon and could ultimately help to 



Time 


Fig. A3. Comparison between the fraction of cloud free, high quality pixels, acquired from MAIAC (black) and MCD43A4 (red). MAIAC observations have been aggregated to 16 day 
periods to meet the temporal scales of MCD43A4. 
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reduce uncertainties of terrestrial vegetation products globally. It 
should be noted that differences between MAIAC and MODIS stan- 
dard processing will likely be considerably smaller in temperate eco- 
systems, as cloud masking is most challenging in tropical regions with 
warm, low clouds that are difficult to detect using conventional tech- 
niques. More research will be required to determine the magnitude of 
this effect in other regions of the globe. Nonetheless, the quality of 
MAIAC demonstrated in tropical ecosystem clearly demonstrates its 
potential for accurate cloud detection and aerosol correction. 

As cloud cover is seasonally dependent, it seems plausible that 
previously reported increases in Amazon greenness during a drought 
could, at least in part, result from reduced cloudiness, rather than 
from a change in vegetation status. More research will be required 
to investigate this possibility but considerable uncertainties in stan- 
dard products found in this study highlight the need for alternative 
processing techniques over tropical regions (Zelazowski et al., 
2011). The results of this study suggest that MAIAC could provide 
an alternative for reassessment of biophysical changes over tropical 
rainforests. 
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Appendix A. Comparison of BRDF normalized products 

In order to investigate the effect of BRDF on MODIS derived NDVI, we 
compared directionally normalized NDVI from the MODIS standard 
product and MAIAC for the same two MODIS tiles (hllv09 and 
hllv08) and time period (2002-2010) (Fig. 1). For that, we used 
MAIAC normalized reflectance BRF n and MODIS NBAR (nadir 
BRDF-adjusted reflectance) product. The BRF n is normalized to the 
fixed geometry of solar zenith angle (0 S ) of 45° and nadir view using 
the RTLS model: 

BRF x RTLS(fl s , fl v , Atp) 

n RTLS(0 s = 45°, 6 V = 0°, A<p = 0°) ’ 

where 9 V and Acp are view zenith and relative azimuth angle, 
respectively. 

The MCD43A4 product is a 16 day composite product which pro- 
vides level-3 gridded surface reflectance at a 0.5-kilometer resolution 
adjusted for view angle effects using the bidirectional reflectance dis- 
tribution function (BRF). MCD43A4 data and corresponding quality 
flags (MCD43A2) were obtained from the EOS data gateway and qual- 
ity filters were applied as stated in Table 1. While MCD43A4 and 
MAIAC BRF n are both BRDF normalized datasets, there are two differ- 
ences between them: First, MCD43A4 is normalized in the view angle 
only, while MAIAC BRF n is normalized in both view and solar zenith 
angles. Second, MCD43A4 is a 16-day composite product based on 
both Terra and Aqua MODIS data, while MAIAC BRF n is daily reflec- 
tance, here from the Aqua platform only. Fig. Al shows NDVI time 
series from the same 50 x 50 km 2 subsets as used before but compar- 
ing MCD43A4 to MAIAC BRF n at monthly time steps. As for the 
non-normalized data, we computed monthly mean values for each 
pixel by averaging all valid (that is cloud screened, and quality 
filtered) pixels acquired during that month and obtained mean and 
standard deviation by averaging these monthly values across each 
50 x 50 km 2 block. Differences observed between the 16 day product 
MCD43A4 and daily MAIAC BRF n were moderate for the most part, 
however, cloud cover caused residual outliers in the 16 day product, 
with some areas more affected than others. As opposed to conven- 
tional composite, MCD43A4 reflectance relies on daily MOD09 obser- 
vations to retrieve BRDF and is also affected by undetected clouds, 


especially if the total number of observations is low. Overall the per- 
formance of the 16 day NBAR composite product was similar to that 
of the 16 day VI product (MYD13). 

Fig. A2 shows NDVI estimates (along z-axis) as a function of time 
during 2002-2010 (x-axis). The y-axis shows the spatial variability of 
NDVI when averaging over different areas, from 2 x 2 to 50 x 50 km 2 . 
The color represents the standard deviation (a) of the spatially aver- 
aged NDVI estimates, the areas correspond to the non-normalized 
data presented in Figs. 8 and 9. MAIAC BRF n observations are shown 
in Figures A (northern, seasonal region), and Figure C (southern, 
non- seasonal part), the corresponding MCD43A4 observations are 
given in Figures B and D. Compared to the non-normalized MAIAC 
product, MAIAC derived BRF n observations showed similar levels of 
uncertainty indicating only a limited effect of BRDF on NDVI reflec- 
tance (compare Figs. 8A and 9A). The MODIS NBAR product yielded 
comparatively overall larger uncertainties across the examined scales 
with individual dates showing considerable offsets, which were likely 
driven by undetected clouds. 

The number of observations shown for the NBAR product was con- 
siderably lower than for MAIAC observations. Fig. A3 shows the pro- 
portion of cloud free pixels from MCD43A4. To match the temporal 
scale of MCD43A4, MAIAC data were aggregated over 16 days. 

In summary, this analysis showed little difference between MAIAC 
BRDF normalized and non-normalized NDVI which agrees with re- 
sults from Lyapustin et al. (this issue). MAIAC NDVI obtained from 
BRF and from BRF n were found to be very similar with only slight re- 
duction in variability. As NDVI is typically high in tropical forests, the 
red reflectance constitutes only a small fraction of the NIR signal, 
thereby canceling most of the anisotropy in NDVI (but not so in the 
reflectance). In addition, the anisotropy gets reduced as the surface 
becomes brighter due to multiple scattering and consequently the 
relative anisotropy of the NIR reflectance is lower than that of the 
red reflectance. 

The results presented in Figs. Al and A2 show that despite addi- 
tional BRDF-based cloud masking in MCD43A4, cloud leak remains 
problematic in this product. While the BRDF retrieval in MCD43 
helped filter residual clouds, when enough clear sky observations 
for reliable BRDF inversion were available, this method does not 
allow the cloud mask error to be separated from the total error bud- 
get. In cases of high cloudiness, residual clouds can lead to wrong 
BRDF retrievals, and, as a result, the MCD43A4 BRDF-normalization 
may increase the total error of surface reflectance (or NDVI) rather 
than decreasing it. It should also be noted that the results presented 
in Fig. Al show monthly-mean values over large areas (2500 km 2 ). 
Acquisitions at finer temporal or spatial scales using MCD43A4 
would be problematic due to a lack of sufficient observations (see 
for instance Figs. A2, A3). 

Appendix B. Supplementary data 

Supplementary data to this article can be found online at http:// 
dx.doi.org/1 0.1 01 6/j.rse.201 2.08.035. 
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